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ABSTRACT 


A  computational  study  of  the  turbulence  model 
parameters  and  their  effect  on  the  flow  field  of  a 
circulation  control  airfoil  was  conducted.  A  Beam-Harming 
algorithm  was  used  to  solve  the  Navier-Stokes  equations  for 
flow  around  a  circulation  control  airfoil.  The  turbulence 
model  was  a  modified  Baldwin-Lomax  model.  The  modification 
was  from  Bradshaw  and  based  on  the  extra  strain  rates 
produced  by  streamline  curvature.  Included  in  the 
modification  is  an  empirical  curvature  correction  constant, 
the  Bradshaw  constant.  The  effect  of  this  empirical 
curvature  correction  on  the  flow  about  an  airfoil  was 
analyzed. 

Using  the  simple  geometry  of  a  cylinder,  the  effects  of 
the  Bradshaw  constant  on  the  flow  field  were  to  be 
monitored.  The  effects  of  the  Bradshaw  constant  on  the  flow 
field  about  the  cylinder  were  then  to  be  compared  to  the 
effects  of  the  Bradshaw  constant  on  the  flow  field  about  a 
circulation  control  airfoil.  By  comparing  these  two 
different  but  yet  similar  flow  fields,  the  application  of 
the  Bradshaw  modification  to  various  geometries  was  to  be 
analyzed.  Comparison  of  the  effects  of  the  Bradshaw 
constant  on  the  two  flow  fields  was  to  provide  an  assessment 
of  the  sensitivity  of  the  Bradshaw  correction  to  different 
situations.  The  questionable  convergence  or  non-convergence 
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of  the  solution  for  the  cylinder  precluded  this  cross 
correlation.  However,  a  more  indepth  investigation,  than 
done  previously,  of  the  effects  of  the  Bradshaw  constant  on 
the  solution  about  an  airfoil  were  analyzed. 

As  the  Bradshaw  constant  is  increased,  the  flow  field 
and  the  flow  coefficients  are  affected.  The  flow  field  in 
the  rear  stagnation  region  develops  from  a  flow  without  any 
vortices  to  a  flow  with  a  primary  and  secondary  vortex  as 
the  Bradshaw  constant  is  increased.  This  twin  vortex  flow 
field,  obtained  at  higher  values  of  the  Bradshaw  constant, 
is  in  close  agreement  with  experimental  results.  Once  a 
high  enough  value  is  reached  for  the  Bradshaw  constant,  the 
flow  coefficients  are  relatively  unaffected  by  any  further 
increase  in  the  curvature  correction.  Por  higher  blowing 
coefficients,  these  same  trends  are  observed. 

When  the  jet  ratios  or  jet  momentum  were  matched  with 
experimental  results,  the  moment  coefficient  provided  good 
results.  However,  the  lift  and  drag  values  were  not  in 
agreement  with  experimental  results.  As  the  accuracy  of  the 
lift  and  drag  values  was  increased,  the  accuracy  of  the 
moment  coefficient  was  degraded. 

This  report  provides  further  insight  into  the  effects 
of  the  Bradshaw-modified  Baldwin-Lomax  turbulence  model. 
However,  the  results  obtained  indicate  that  a  better 
turbulence  model  must  be  employed  before  CFD  can  be  used  to 
accurately  predict  the  flow  around  a  circulation  control 
airfoil . 
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COMPUTATIONAL  INVESTIGATION  OF  CIRCULATION 


CONTROL  TURBULENCE  MODELING 

I.  INTRODUCTION 

How  can  we  increase  the  airfoil's  lift?  What  can  be 
done  to  maximize  the  lift  coefficient?  Over  the  years, 
researchers  have  been  continiously  investigating  these 
questions.  Numerous  concepts  have  been  introduced  to  answer 
these  questions.  One  concept  generating  a  lot  of  interest 
in  the  last  few  years  is  the  circulation  control  airfoil. 

Circulation  control  augments  an  airfoil's  lift  by 
blowing  a  thin  jet  over  a  rounded  trailing  edge.  By  taking 
advantage  of  the  Coanda  effect,  the  jet  delays  separation, 
forces  the  stagnation  point  to  the  lower  surface,  and 
increases  circulation  around  the  airfoil. 

The  Coanda  effect  is  the  ability  of  a  tangential  jet  to 
stay  attached  to  a  curved  surface  through  large  turning 
angles  (6:53).  The  entrainment  of  the  flow  between  the  wall 
and  the  jet  creates  a  low  pressure  region  causing  the  jet  to 
attach  to  the  surface.  The  suction  generated  by  the 
entrainment  is  sufficient  to  overcome  the  centrifugal 
forces,  thus  causing  the  jet  to  bend  with  the  wall.  The  jet 
finally  separates  from  the  wall  due  to  viscosity  induced 
momentum  losses  (6:54). 
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The  main  difference  between  the  circulation  control 
airfoil  and  a  conventional  airfoil  is  the  trailing  edge 


configuration.  Instead  of  a  sharp  trailing  edge,  the 
circulation  control  airfoil  employs  a  rounded  trailing  edge. 
Thus,  the  Kutta  condition,  imposed  on  airfoils  with  sharp 
trailing  edges,  is  not  imposed  on  a  circulation  control 
airfoil  (32:1).  Instead,  the  location  of  the  stagnation 
point  is  a  function  of  the  momentum  of  the  jet  and  the 
freestream  flow  conditions.  Thus,  by  varying  the  jet 
momentum  the  airfoil's  lift  can  be  varied  (17:1).  Figure  1 
displays  the  geometry  of  a  typical  circulation  control 
airfoil . 


With  this  capability,  a  circulation  control  airfoil  can 
obtain  higher  lift  coefficients  and  lower  aircraft  stall 
speeds  than  conventional  airfoils.  This  allows  several 
benefits  from  reduced  runway  lengths  to  increased  payload 
weights  (32:1).  However,  the  circulation  control  airfoil 
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experiences  more  drag  than  a  typical  airfoil.  The  rounded 
trailing  edge  creates  a  larger  wake  region  than  a  sharp 
trailing  edge  airfoil.  The  wake  region  can  be  reduced  by 
using  a  flap  type  system  to  create  a  sharp  trailing  edge 
when  the  additional  lift  is  not  needed.  The  added  weight 
and  complexity  of  pressurized  ducts  is  another  disadvantage 
of  the  circulation  control  airfoil  concept  (31:3). 

The  circulation  control  concept  has  many  applications. 
Short  takeoff  and  landing  (STOL)  aircraft  can  benefit  by 
reducing  forward  speed  while  increasing  lift.  Circulation 
control  airfoils  have  been  investigated  for  application  on 
the  X-wing  stopped  rotor  concept  aircraft.  The  no-tail 
rotor  (NOTAR)  helicopter  design  uses  circulation  control  to 
produce  the  anti-torque  forces  required  for  stable  flight. 
These  are  just  a  few  present  applications  of  the  circulation 
control  airfoil  concept  (17:16). 

1 • 1  Background 

About  1910,  Henry  Coanda  was  trying  to  deflect  the 
engine's  exhaust  away  from  his  aircraft's  wooden  wings. 
Instead  of  deflecting  the  exhaust,  the  exhaust  adhered  to 
the  deflection  plates.  This  led  to  two  important 
discoveries  by  Henry  Coanda:  1)  the  "Coanda  effect";  and 
2)  burning  wings  don't  produce  much  lift  (32:3). 

Since  then  a  tremendous  amount  of  time  has  been  spent 
on  researching  the  circulation  control  airfoil  concept.  The 
National  Gas  Turbine  Establishment  in  England  started  doing 
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extensive  studies  into  the  benefits  of  circulation  control 
airfoils  about  1960  (7).  About  the  same  time,  V.E.  Lockwood 
did  some  experimental  work  with  a  blowing  slot(s)  on  a 
cylinder  (15).  Throughout  the  1970s  and  up  to  present  day, 
the  David  Taylor  Naval  Ship  Research  and  Development  Center 
(DTNSRDC)  has  compiled  an  extensive  database  on  circulation 
control  airfoils  (9).  In  1986,  NASA  held  a  circulation 
control  workshop  at  Ames  Research  Center  (18).  This 
conference  presented  the  latest  developments  in  experimental 
and  computational  research  and  potential  applications  for 
circulation  control  airfoils.  A  good  overview  of  the 
conference  was  compiled  by  Nielsen  and  Biggers  (17). 

As  Nielsen  and  Biggers  state,  the  computational  codes 
"are  promising,  but  as  yet  inadequate  for  a  priori  modeling 
of  the  Coanda  flows."  This  inadequacy  is  due  to  the  lack  of 
an  accurate  model  for  the  flow  around  the  jet  slot  and  the 
Coanda  surface.  How  does  the  jet  blowing  interact  with  the 
freestream  boundary  layer?  What  complications  does  the 
rounded  Coanda  surface  impose  on  the  boundary  layer?  And, 
what  type  of  corrections  must  be  implemented  to  account  for 
these  effects?  Even  though  researchers  have  been  able  to 
get  good  results,  they  agree  that  these  questions  must  be 
answered  before  the  computational  algorithms  can  be  used  to 
predict  circulation  control  performance. 

Pulliam  et  al  (20)  used  the  thin-layer  Navier-Stokes 
equations  in  generalized  coordinates  to  compute  the  flow 
around  a  circulation  control  airfoil.  They  used  an  implicit 
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approximate  factorization  technique  with  the  zero  equation 
Baldwin-Lomax  turbulence  model  to  solve  the  governing 
equations.  In  addition,  Pulliam  et  al  used  a  correction  to 
the  Baldwin-Lomax  model  suggested  by  Baldwin  to  account  for 
the  curvature  effects  of  the  rounded  trailing  edge.  A 
spiral  grid  was  used  to  compute,  instead  of  model,  the  flow 
at  the  slot  exit.  Even  though  the  results  were  comparable 
to  experimental  results,  the  curvature  correction  had  to  be 
adjusted  to  "match"  the  computational  results  with  the 
experimental  results.  As  long  as  this  "matching"  process  is 
required,  the  code  can  not  be  used  to  predict  the 
performance  of  circulation  control  airfoils. 

To  account  for  the  effects  of  curvature,  Shrewsbury 
(26)  and  Williams  (31)  modified  the  Baldwin-Lomax  turbulence 
model  with  an  empirical  curvature  correction  factor  put 
forth  by  Bradshaw  (6:68-71).  Both  researchers  used  this 
modified  turbulence  model  in  conjunction  with  Beam-Warming 
related  algorithms  to  model  the  flow  around  a  circulation 
control  airfoil.  Again,  both  researchers  obtained 
reasonable  results,  but  only  after  adjusting  the  curvature 
correction  to  "match"  the  computational  results  with  the 
experimental  results. 

As  Nielsen  and  Biggers  stated,  these  computational 
codes  do  have  promise.  The  Navier-Stokes  solvers  used  by 
Pulliam,  Shrewsbury,  and  Williams  do  provide  reasonable 
results.  However,  the  trial  and  error  process  associated 
with  "matching"  the  computational  results  with  the 
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experimental  results  through  the  varying  of  the  empirical 
turbulence  corrections  must  be  eliminated.  If  computational 
methods  are  to  be  used  in  the  research  and  design  of 
circulation  control  airfoils,  an  approach  must  be  developed 
that  does  not  require  "matching”.  Two  courses  of  action  are 
needed  to  eliminate  this  trial  and  error  process: 

1)  experimental  efforts  need  to  concentrate  on  analyzing  all 
the  aspects,  especially  the  turbulent  makeup,  of  the  complex 
flow  field  caused  by  the  jet  and  Coanda  surface; 

2)  computational  programmers  need  to  develop  a  broader 
knowledge  about  the  effects  of  the  turbulence  correction 
factors.  With  this  increased  knowledge,  the  need  to  "match" 
computational  results  with  experimental  results  may  be 
eliminated.  Thus,  the  computational  codes  can  then  be  used 
to  predict  the  performance  of  circulation  control  airfoils. 

1 • 2  Objective 

The  purpose  of  this  study  is  to  further  the  efforts 
associated  with  the  second  course  of  action  described  above. 
Instead  of  validating  the  capability  of  the  code  to  model  a 
circulation  control  airfoil,  as  has  been  done  previously, 
this  research  examines  the  correlation  between  the 
computational  results  and  the  Bradshaw  curvature  correction. 
Using  the  Beam-Warming  algorithm  and  modified  Baldwin-Lomax 
turbulence  model  employed  by  Williams,  a  systematic  approach 
to  the  use  of  the  curvature  correction  factor  is  examined. 
Instead  of  just  "matching"  computational  results  with 
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experimental  results,  this  study  will  investigate  the  effect 
of  the  Bradshaw  curvature  correction  on  the  computed  flow 
field.  In  addition  to  examining  the  variation  in  the 
coefficients,  how  does  the  actual  flow  field  change  as  the 
Bradshaw  correction  is  increased?  Can  any  trends  in  the 
effect  of  the  Bradshaw  correction  be  revealed  that  will 
allow  the  computational  codes  to  be  used  to  predict 
circulation  control  flow  fields?  To  generalize  this 
correlation  for  a  variety  of  flow  geometries,  the  algorithm 
is  applied  to  the  "simple"  geometry  of  a  two-dimensional 
circular  cylinder.  The  effect  of  the  Bradshaw  correction  on 
the  flow  field  about  the  circular  cylinder  will  then  be 
compared  to  the  effect  of  the  Bradshaw  correction  on  the 
flow  field  about  the  airfoil.  By  comparing  the  trends  for 
both  flow  fields,  an  assessment  of  the  curvature  corrections 
effect  for  various  flow  geometries  can  be  ascertained. 
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II.  MIALYSIS 


The  complex  nature  of  the  flow  around  a  circulation 
control  airfoil  requires  the  use  of  the  Navier-Stokes 
equations.  As  a  group,  the  Navier-Stokes  equations  comprise 
the  conservation  of  mass  (continuity),  conservation  of 
momentum  (Newton's  second  law),  and  conservation  of  energy 
(first  law  of  thermodynamics).  The  following  subsection 
will  present  the  derivation  of  the  strong  conservation 
vector  form  of  the  Navier-Stokes  governing  equations.  The 
remaining  subsections  discuss  the  numerical  implementation 
of  the  Navier-Stokes  equations,  the  application  of  the 
boundary  conditions,  the  modifications  made  to  the  Baldwin- 
Lomax  turbulence  model,  and  the  generation  of  the  grids. 

2 . 1  Governing  Equations 

For  unsteady,  compressible,  viscous  flows  the  governing 
equations  are  (19:11) 

Conservation  of  Mass  - 

-|£+?*<p?)=0  (1) 

Conservation  of  Momentum  - 
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Conservation  of  Energy  - 


+v-Et?=  |e + -?•<? 


The  total  energy,  El  ,  is  expressed  as 


Ec  =  p(e  +  -iv'2) 


For  the  case  of  a  Newtonian  fluid  the  stress  tensor,  o,  is 

n=-px+i  (5 

where  the  viscous  stress  terms  are  defined  as 

&  ,  dUi  duj.  duk 

The  viscosity  functions  are  calculated  using  Stokes’ 


hypothesis  (23:60) 


and  Sutherland's  formula  (23:328) 


r3/2 

r+Cj 


Closure  of  the  Navier-Stokes  equations  is  provided  by  the 
perfect  gas  law 


p  =  pf?r=pe(Y-D 
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and  Fourier's  law  of  heat  conduction 


<f=  -kVT  <10> 

With  the  approximation  of  constant  Prandtl  number,  Pr,  the 
coefficient  of  thermal  conductivity,  k,  is 

*  =  (11) 

Pi 

For  a  two-dimensional  rectilinear  system,  the 
Navier-Stokes  equations  become  (16:5) 

Conservation  of  Mass  - 

dp  d(pu)  d(py)  _0  (12) 

5t  dx  dy 

Conservation  of  Momentum  - 

d(pu)  ,  d(pu2*p-xxx)  dipuv-x^)  _n  (13) 

at  dx  dy 

3(pv)  |  djpuv-x^  _  dtp (14) 

3t  dx  dy 

Conservation  of  Energy  - 

dEt  t  d(Ecu+pu-uxxx-vx 

dt  a,-  dx  ,  (15) 

d(Etv+pv-ux  „-vx  „*qy) 

Note,  the  body  forces  are  assumed  negligible  and  the  heat 
generation  term  is  assumed  constant  with  respect  to  time. 
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To  facilitate  the  application  of  these  equations  to  a 
general  curvilinear  coordinate  system,  the  physical  domain 
(x,y)  must  be  transformed  to  the  computational  domain  ((,?)). 
The  transformation  relating  the  physical  and  computational 
domain  is 


Z*Z(x,y)  (16) 

r\  =11  < x,y ) 

Through  application  of  the  chain  rule,  the  physical  domain 
derivatives  become 


d{  )  _  35  3(  )  t  dr\  3(  ) 

dx  dx  3$  &x  3ii 


(17) 


d(  )  .  35  3(  )  3r)  3(  ) 

dy  dy  35  3y  3i) 


(18) 


After  applying  the  chain  rule  and  non-dimensional izing  (see 
Appendix  A),  the  Navier-Stokes  equations  in  matrix  form  for 
a  general  curvilinear  coordinate  system  are 


dU  P  dF  +  t  dG 
dt  ^x35  Qy35 


=  0 


(19) 


where 


U  = 


P 

P  U 

pv 

\Et/ 


(20) 
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(21) 


P  u 

P  ua-t„ 

PUV-Xxy 

(Et+p)  u-irc^-vr 


pv 

puv-x^ 

Pv2~xyy 

(Et+p)  v-irz  -vz 


(22) 


The  first  equation  is  the  conservation  of  mass;  the  center 
equations  are  the  x  and  y  momentum  equations  respectively 
and  the  final  equation  is  the  conservation  of  energy. 

The  turbulence  is  included  through  the  turbulent  eddy 
viscosity,  e,  and  the  turbulent  coefficient  of  thermal 
conductivity,  kt.  With  the  inclusion  of  the  turbulence, 
Stokes'  hypothesis  becomes 

<23) 

and  the  turbulent  coefficient  of  thermal  conductivity  is 
obtained  by  defining  the  turbulent  Prandtl  number,  Prt, 


(24) 


Sutherland's  formula,  Eq  (8),  is  used  to  calculate  the 
molecular  viscosity,  p;  the  turbulent  eddy  viscosity  is 
calculated  with  the  turbulence  model  (sec  2.4);  the 
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turbulent  Prandtl  number  is  assumed  constant.  With  these 
definitions,  the  viscous  stress  and  heat  conduction  terms 
become  (3:86) 

****-3*tU*+*t<U*+Vy)  <2S> 


Tw=-|ic(Uy+Vx) 

(26) 

t =  -3A,  tvy* k  t(ux*vy) 

(27) 

**=-<*+*c)|£ 

(28) 

qy  =  -  (,k+kt) 

(29) 

Using  the  definition  of  the  transformation  Jacobian,  J 

T-  -F  n  -F  n 

duly)  (yn* 

(30) 

the  Navier-Stokes 

conservation  form 

equations  are  rewritten  in  strong 

(29:7) 

where 

3ff+  df  d<2. 

(31) 

U 

(32) 
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(33) 


P^-jaxF  +  lyG) 

<?=-i(t|  *F-M)yG)  (34) 

By  separating  the  viscous  and  non-viscous  terms,  Eq  31  is 
rewritten  as 


dU  d ^  _  dV  dW  (35) 

dt  dl  dr\  ~  dt  dr) 

This  revision  separates  the  functional  relationships  since 

F^f(U)  and  61=f(U)  <36> 

but 

V‘f(8,Gv8j  and  W*f(8,8(,u^)  (”) 

A  further  refinement  separates  the  derivative  relationships 

SO  dF,(0)  60,(0) 

Jt+  3{  +  3n  (38) 

dv2(0,0,)  dv2(0,0.)  6W,(U,0.)  dW2(0,0.) 

— M — +  H  * — Si — +  6n 

Note  V=Vj+Vj  and  WsW^+Wj.  By  substituting  Eqs  (25) -(27)  in 
for  the  shear  stresses  and  using  Stokes'  hypothesis,  Eq 
(23),  the  matrices  become 
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1 

J 


p 

pu 

pv 

\Et> 


ft  _  1 


puc  1 
PUL2C  +  ^XP 
pvue+£yp 
^  <£t+p)  UCJ 

P^c 

PUVC+T\XP 
PWc+V  yP 
{  (Et+p)  vc) 


(39) 


(40) 


(41) 


*1  = 


_1 

J 


b1u^b2vl 

^u^b^ 

[b^u^ + b2  ( vu{  +  uv{ )  +  i>3  wj  +  jb4  Tj 


(42) 


v,  =  -4 


1 

J 


c3i^  +  cAvn 
[c^  +  c2uv^  +  c3  vu„  +  C.w,,  +  C5rj 


(43) 


"i= 


_1 

J 


CiUc  +  C3Vj 

c2u^cAv% 

^ctuu?  +  c2vu?  +  C3liV{  +  c4w(  +  c5rtJ 


(44) 
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1 

1 

1 

r  0  > 

■ 

1 

II 

<4‘Vd3v<t 

(45) 

i 

^ui^  +  djlvi^  +  uv,)  +d3w,  +  d4T,J 

1 

where  the  contravariant  velocities  are 

1 

Ue"«*u+?yv’ 

(46) 

1 

Vc  *  T)  XU  +  T|  yV 

(47) 

1 

and  the  viscous 

coefficients  are 

1 

£>1=  <|»+e)  5  *+{’,) 

(48) 

1 

1 

bZ  =  ±(V*t)l£y 

(49) 

1 

^3 =  (§i  +  e)  ( 5  *  +  5  y) 

(50) 

1 

■ 

=c’{~pI*-pFt) 

(51) 

1 

1 

C^-iVL  +  e)  (-l-^x^X+^y^y) 

(52) 

1 

1 

1 

c,  =  -(|t  +  e)  (-ICxVM*) 

(53) 
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(54) 


c3=  (n  +  e)  (C,ny--|5y'i*) 


c4  =  -< in-t)  <C*n*+-§-Cy»iy> 

(55) 

(56) 

c?x=  <H+e) 

(57) 

d2  =  -|(H+e)i)Jtf)y 

(58) 

d3=  (Ji+E)  <flx+-jT)2y) 

(59) 

d*sC>{~pZ  +  ’W~)  {1'2*+v'2y) 

(60) 

2 . 2  Numerical  Implementation 

The  numerical  technique  used  to  solve  the  Navier-Stokes 
equations  is  based  on  the  Beam  and  Harming  approximate 
factorization  algorithm  (5,30)  with  a  Baldwin-Lomax 
Turbulence  Model  (4).  The  code  was  developed  by  Dr.  Miguel 
R.  Visbal  of  the  Air  Porce  Wright  Aeronautical  Laboratories 
(28,29).  Steve  Williams  added  some  modifications  to  the 
code  to  adapt  it  to  circulation  control  airfoils  (31). 

The  algorithm  solves  the  strong  conservation  form  of 
the  Navier-Stokes  equations  (Eq  38).  The  solution  algorithm 
for  first  order  Euler  time  differencing  is  written  as 


17 


a«  a«j 


3*1  dtl2 


-At [4-  (^i-Vi-Va)  “♦-g-  ((%-Nk-l^)  fl] 


^v"i  ^  ^  ’“Sf 

The  Jacobian  matrices  A,  B,  R,  and  S  are  (see  Appendix  B) 


(SL) 


A  = 


i4 

so 

ea 

avi 

s.dH? 

dfy 

(62) 


(63) 


and  n  is  the  temporal  index.  The  equation  is  solved  in 
three  steps: 

1)  along  lines  of  constant  ^ 


[  I  +  A  t  ( ^  °  ]  D  D= 
dl  dt2 


-A t[-^  *  +  (^Wx-W2)  D] 

2)  along  lines  of  constant  £ 

[x  +  At(^^--^^)]  Ai7"=r>» 

«i  9n 


3)  update  solution 


(64) 


(65) 


=  ffn  +  Atfn 


(66) 


The  spatial  derivatives  are  computed  using  second  order 
central  differences.  This  setup  requires  the  solution  of  a 
block  tridiagonal  system  for  steps  1)  and  2)  along  each 
coordinate  line,  where  each  block  has  dimensions  4x4. 
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To  damp  any  numerical  instabilities,  artificial 
viscosity  terms  are  added  to  the  algorithm.  An  explicit 
fourth  order  term  is  added  to  the  right-hand  side  of  Eq  (61) 

D.=  -±a,&t(l*t  +  6\)Un  <”> 

In  addition,  second  order  implicit  damping  terms  are  added 
to  enhance  stability.  The  implicit  terms  are  inserted 
within  the  applicable  operators  on  the  left-hand  side  of 
Eq  (61) 

JT  <68> 

JT  (69> 

The  6  operators  are  central  difference  operators  in  the 
respective  directions  and  of  the  respective  order. 

Local  time  stopping  is  introduced  to  accelerate 
convergence  to  a  steady  state  solution.  The  time  step  is 
computed  by 

AC-CA  (70) 

where  C  is  the  Courant  number.  Appendix  C  contains  a  list 
and  description  of  the  required  program  inputs. 

Convergence  of  the  code  was  determined  by  monitoring 
the  oscillations  of  the  lift  and  drag  coefficients.  A 
solution  was  assumed  valid  when  the  amplitude  of  the 
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oscillations  was  less  than  0.1  percent  for  both  the  lift  and 
drag  coefficients.  The  convergence  criteria  was  modified 
for  the  cylinder  geometry  at  high  Reynolds  numbers  since  a 
steady-state  solution  was  not  expected  due  to  the  presence 
of  vortex  shedding.  The  convergence  criteria  for  the 
cylinder  at  high  Reynolds  numbers  was  to  allow  only  a  one 
percent  scatter  about  a  periodic/quasi -periodic  pattern. 

2 . 3  Boundary  Conditions 

Boundary  conditions  must  be  specified  at  the  grid 
exterior  boundary,  the  grid  cut  line,  the  airfoil  surface, 
and  the  jet  slot.  The  exterior  boundary  of  the  grid  is 
divided  into  an  inflow  and  outflow  portion.  Typical 
freestream  conditions  are  used  along  the  inflow 

u*{7„cos(a)  (71) 

v*t7„sin(a)  (7*) 


P-P- 

P-P. 

Along  the  outflow,  subsonic  velocity  is  assumed  with 


(73) 

(74) 


dv 

"Si 


®  o 


(75) 


(76) 
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p-p. 


(77) 


|H-0  (78) 

dx 

The  grid  overlapped  two  points  on  both  sides  of  the  cut 
line.  Along  these  overlapped  lines  periodic  boundary 
conditions  are  enforced. 

On  the  airfoil  surface,  except  the  jet  slot,  no-slip 
and  adiabatic  wall  boundary  conditions  are  imposed: 

U—  V=  0  (79) 

-gH  =0  (80) 


dT 

IE 


=  0 


(81) 


The  jet  slot  is  assumed  convergent  with  total 
temperature  and  pressure  constant  across  the  exit.  The  flow 
out  the  slot  is  assumed  to  be  isentropic.  The  pressure  at 
the  slot  is  extrapolated  with  the  use  of  Eq  (82) 


(82) 


Using  the  jet  static  pressure  calculated  from  Eq  (82)  and 
the  input  variable  of  jet  total  pressure  over  freest ream 
static  pressure  (see  Appendix  C),  the  ratio  of  jet  static 
pressure  to  jet  freestream  pressure  can  be  determined. 
Using  the  isentropic  relationship. 
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(83) 


the  jet  Mach  number  can  be  calculated.  With  the  jet  Mach 
number  and  the  ratio  of  jet  total  temperature  to  freestream 
static  temperature  (input  variable),  the  jet  static 
temperature  can  be  calculated  with  the  iaentropic 
relationship 


(84) 


With  the  jet  Mach  number  and  static  temperature,  the  jet 
velocity  components  can  be  calculated 

U^M(yRT)  1/2  COS<|>  (85) 

v=M(yRT)  1/2sin$  (86) 

where  4  is  the  angle  of  the  jet  midline  with  the  coordinate 
axes.  Note  the  flow  at  the  exit  is  limited  to  sonic 
velocity  due  to  chocking  (31). 

The  input  values  of  jet  pressure  and  temperature  ratio 
and  the  slot  height  dictate  the  jet  momentum  coefficient. 
The  jet  momentum  coefficient  is  defined  as  (31:17) 


E1Z1 

Q.c 


(87) 


Note,  the  jet  velocity  is  determined  assuming  isentropic 
expansion  out  the  jet  slot. 
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The  drag  calculation  is  modified  due  the  presence  of 
the  jet  thrust.  The  modified  drag  coefficient  is  calculated 
by  subtracting  the  momentum  coefficient  from  the  friction 
and  pressure  drag  (31:18) 

cd  <88) 

2 . 4  Turbulence  Model  and  Point  of  Transition 

The  accuracy  of  the  computational  solution  for  flow 
around  a  circulation  control  airfoil  is  heavily  dependent  on 
correctly  predicting  the  point  of  boundary  layer  detachment 
from  the  Coanda  surface.  The  turbulence  model  strongly 
influences  the  location  of  the  separation  point.  Thus, 
without  a  valid  turbulence  model  the  computational  solution 
is  not  accurate. 

As  discussed  previously,  the  turbulence  model  most 
widely  used  for  circulation  control  modeling  is  the  Baldwin- 
Lomax  model  (4).  This  zero  equation  model  is  widely  used 
due  to  its  simplicity  and  versatility.  But,  this  model 
requires  modifications  to  handle  the  flow  around  the  Coanda 
surface. 

The  rounded  Coanda  surface  can  have  a  surprisingly 
large  effect  on  the  turbulent  boundary  layer.  Streamline 
curvature,  caused  by  the  Coanda  surface,  is  a  dominant 
source  of  additional  strain  rates  in  turbulent  boundary 
layers.  The  extra  rates  of  strain  can  have  both  a 
stabilizing  and  destabilizing  effect  on  the  flow.  Balancing 
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the  centrifugal  and  pressure  forces  on  a  fluid  element 
results  in  a  stabilizing  effect  in  regions  where  the 
velocity  increases  normal  to  the  surface,  3U/3n  >0.  In 
regions  where  the  velocity  decreases  normal  to  the  surface, 
3U/3n  <  0,  streamline  curvature  has  a  destabilizing  effect 
on  the  flow.  These  instabilities  can  increase  or  decrease 
the  turbulent  transport  of  the  flow.  The  streamline 
curvature  of  a  flow  field  can  be  an  order  of  magnitude  more 
important  to  turbulent  boundary  layer  development  than  the 
normal  pressure  gradient  (6). 

To  account  for  this  behavior,  Bradshaw  (6)  recommends 
adjusting  the  mixing  length.  The  Baldwin-Lomax  model 
calculates  the  turbulent  eddy  viscosity  of  the  inner  layer 


(|*t>  lnn«“ P“J2 

To  account  for  the  streamline  curvature,  the  length  is 
multiplied  by  a  curvature  correction  factor,  F, 


F  =  1-  0S  (90> 

The  parameter  S  represents  the  ratio  of  extra  strain  rates 
induced  by  the  streamline  curvature  to  the  inherent  strain 
rates 


(91) 


8  is  an  empirical  curvature  correction  constant  which  varies 
from  case  to  case  but  is  of  the  order  of  ten. 
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Adverse  pressure  gradients  normal  to  the  flow  have  a 
big  impact  on  boundary  layer  transition.  Laminar  flows  can 
not  withstand  any  appreciable  adverse  pressure  gradients. 

The  adverse  pressure  gradient  tends  to  amplify  any  initial 
velocity  oscillations,  thus  increasing  the  tendency  to 
transition  to  fully  turbulent  flow.  The  appearance  of  the 
low  frequency  Tol lmien-Schlichting  waves  in  the  laminar 
boundary  layer  is  the  first  indication  of  boundary  layer 
instability.  As  t.ij  velocity  oscillations  grow,  the  laminar 
boundary  layer  quickly  degenerates  into  a  fully  turbulent 
boundary  layer.  Even  though  the  phrase  “point  of 
transition"  is  used  throughout  the  literature,  the 
transformation  from  a  laminar  to  a  turbulent  boundary  layer 
occurs  over  a  finite  distance  (23:445-500). 

This  study  took  a  conservative  approach  towards 
specifying  the  transition  from  laminar  to  turbulent.  Using 
the  Pohlhausen  method  based  on  the  elliptic  airfoil's 
slenderness  ratio  (23:492-499),  the  point  of  instability  was 
determined.  Then  the  pressure  field  between  the  front 
stagnation  point  and  the  point  of  instability  is  examined. 
The  minimum  pressure  point  located  in  this  region  is  then 
specified  as  the  “point  of  transition".  This  criteria  leads 
to  an  increase  in  calculated  drag  when  compared  to 
experimental  results  since  the  actual  “point  of  transition" 
is  downstream  of  the  minimum  pressure  point  and  point  of 
instability,  in  general  (23:489-500). 
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2 . 5  Grid  Generation 

Due  to  the  rounded  trailing  edge  and  variable 
separation  point  location,  circulation  control  airfoils  are 
not  conducive  to  the  use  of  C-grids.  However,  an  O-grid 
conforms  to  the  Coanda  surface  and  the  cut  line  does  not 
have  to  be  located  in  the  wake  region.  For  both  the 
circular  cylinder  and  the  airfoil,  O-grids  were  used  with 
the  cut  line  at  the  leading  edge  of  the  model.  All  grids 
were  generated  with  GridGen-2D,  a  software  package  developed 
by  General  Dynamics  for  the  Air  Force. 

Once  the  surface  contour,  cut  line,  and  exterior 
boundary  were  defined,  Vinikour  stretching  was  used  to 
cluster  the  grid  points  around  the  critical  regions,  i.e. 
jet  slot  and  Coanda  surface.  The  grid  lines  were  generated 
using  a  Soni  TFI  algebraic  solver  with  further  refinement 
accomplished  using  an  elliptic  solver.  GridGen-2D  allows 
several  variations  on  the  type  of  elliptic  solver.  The 
grids  used  in  this  study  were  developed  using  Poisson's 
equation  with  a  Sorenson  weighting  function  and  with 
orthogonality  enforced  at  the  surface.  The  grid  around  the 
entire  airfoil  is  shown  in  Figure  2;  Figure  3  is  a  closeup 
of  the  jet  slot  and  Coanda  surface. 

A  major  headache  in  circulation  control  airfoil  grid 
generation  is  the  discontinuous  surface  geometry  at  the  jet 
slot.  The  sharp  corners  dictate  a  high  degree  of  grid 
skewness.  In  addition  to  the  various  techniques  used  with 
the  GridGen-2D  elliptic  solver,  a  hyperbolic  solver  was 


investigated  to  minimize  the  grid  skewness  (12).  However, 
neither  solver  was  able  to  drastically  reduce  the  skewness, 
see  Figure  4.  The  poor  grid  metrics  caused  by  the  skewness 
can  have  a  degrading  influence  on  the  solution.  A 
qualitative  value  for  this  degradation  is  unknown. 

For  the  one  inch  diameter  circular  cylinder,  the  final 
grid  had  207  points  azimuthal ly  and  151  points  normal  to  the 
surface.  To  insure  accurate  modeling  of  the  boundary  layer, 
normal  spacing  at  the  wall  was  0.00001  chords.  The  far- 
field  boundary  was  circular  with  a  diameter  of  30  chords. 

In  order  to  compare  answers  with  previous  analysis,  the 
103RE  (also  known  as  the  103XW)  airfoil  used  by  Williams 
(31)  was  used  for  this  study.  However,  the  surrounding  grid 
was  modified  slightly.  The  number  of  points  around  the 
airfoil  was  increased  from  176  to  211  to  provide  better 
resolution  over  the  Coanda  surface.  In  comparison  with  the 
flow  field  around  the  cylinder,  the  streamlined  flow  around 
the  airfoil  reduces  the  threat  of  a  reflected  pressure  wave. 
Therefore,  the  far-field  diameter  was  reduced  to  10  chords. 
Eighty  points  were  distributed  between  the  airfoil  surface 
and  the  far-field  boundary  with  an  initial  spacing  at  the 
wall  of  0.000005  chords. 
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Lgure  2.  Inner  Region  of  Airfoil  Grid 


Jet  Exit/Coanda  Surface  Region  of  Airfoil  Grid 


III.  RESPLTS  AND  DISCOSSIOH 


To  investigate  the  complexities  of  flow  around  a 
circulation  control  airfoil,  two  different  geometries  were 
analysed  *  a  circular  cylinder,  and  an  airfoil.  The  flow 
around  both  geometries  was  subjected  to  a  Beam-Harming 
computational  solution  of  the  Navier  Stokes  equations. 

This  chapter  discusses  the  results  obtained  for  both 
geometries.  The  first  portion  of  the  chapter  discusses  the 
results  obtained  with  the  circular  cylinder.  The  results 
for  the  103RE  cross  section  experimental  circulation  control 
helicopter  rotor  follow  the  cylinder  discussion. 

3.1  Circular  Cylinder 

The  flow  around  a  circular  cylinder  is  characteristic 
of  flow  past  a  "bluff  body"  (11:30).  For  a  bluff  body,  the 
flow  separates  well  ahead  of  the  rear  trailing  edge,  thus 
producing  a  large  wake  region.  The  exact  location  of  the 
separation  point  and  the  nature  of  the  flow  in  the  wake 
region  are  highly  dependent  on  the  flow  Reynolds  number. 
Appendix  D  contains  a  breakdown  of  the  flow  regimes  around  a 
circular  cylinder  for  a  range  of  Reynolds  numbers. 

The  multifarious  nature  of  the  flow  around  a  circular 
cylinder  requires  a  systematic  analysis.  The  first  step  was 
to  analyze  the  flow  around  a  cylinder  without  a  jet  slot. 
Using  a  grid  of  155  x  76,  the  code  converged  to  a  steady 
state  solution  for  a  Reynolds  number  of  40  with  a  freestream 
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Mach  number  of  0.2.  A  plot  of  the  velocity  vectors  in  the 
rear  separated  wake  region  is  shown  in  Figure  5.  The 
computea  drag  coefficient  was  1.53,  which  is  in  good 
agreement  with  experimental  values  and  computational  results 
(28).  This  case  served  to  validate  the  code  for  this 
Reynolds  number  regime  and  was  used  as  a  starting  point  for 
cases  with  jet  slot  blowing  and  higher  Reynolds  numbers. 

A  slot  was  then  added  at  the  top  of  the  cylinder,  see  Figure 
6.  The  slot  height  was  chosen  as  0.01  to  match  a  model 
currently  being  used  at  AFIT.  Three  fourths  of  the 
cylinder,  from  the  rear  horizontal  plane  clockwise  to  the 
top  vertical  plane,  had  a  constant  radius  of  0.5.  Due  to 
the  presence  of  the  slot,  the  upper  right  quadrant  of  the 
cylinder  had  an  increasing  radius  of  curvature  from  0.49  at 
the  top  to  0.5  at  the  rear  trailing  edge.  A  155  x  155  grid 
was  generated  with  clustering  around  and  downstream  of  the 
slot  exit. 

The  flow  was  computed  for  a  Reynolds  number  of  40  and  a 
Mach  number  of  0.2  without  blowing.  The  presence  of  the 
slot  discontinuity  caused  several  disturbances  in  the  flow 
field.  Downstream  of  the  slot  several  small  vortices  were 
apparent  along  with  some  alternating  recirculation  regions. 
In  addition,  the  symmetric  counter  rotating  vortices, 
characteristic  of  a  cylinder  without  a  slot,  were  no  longer 
symmetric  for  the  slotted  cylinder.  The  vortices,  along 
with  the  rear  stagnation  streamline,  were  skewed  downward 
and  the  strength  of  the  upper  vortex  increased  while  the 
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slot,  Re=40 


strength  of  the  lower  vortex  decreased.  (Figures  8  and  9 
show  the  same  general  trends  observed  for  this  converged 
Re=40  case).  These  disturbances  resulted  in  an  increase  in 
the  drag  coefficient  from  1.53  without  the  slot  to  1.59  with 
the  slot.  The  addition  of  the  slot  also  caused  an  increase 
in  the  number  of  iterations  required  for  convergence. 

Without  the  slot,  the  solution  converged  in  approximately 
900  iterations.  With  the  slot,  convergence  was  not  achieved 
until  approximately  2400  iterations.  The  initial  condition 
used  for  both  the  slotted  cylinder  and  the  clean  (unslotted) 
cylinder  was  a  uniform  flow  field  throughout  the  grid. 

With  a  uniform  flow  field  throughout,  the  code  was  run 
for  several  hundred  iterations  (*  500)  before  the  boundary 
conditions  at  the  slot  were  changed  to  initiate  blowing. 

The  new  boundary  conditions  corresponded  to  a  momentum 
coefficient  of  approximately  0.01.  Once  the  jet  blowing  was 
initiated  the  solution  tended  to  diverge.  Abnormal  velocity 
vectors  were  being  computed  in  the  vicinity  of  the  slot.  A 
representative  example  of  the  flow  field  around  the  slot  is 
diagramed  in  Figure  7.  The  next  attempt  allowed  the 
non-blowing  case  to  converge  completely  (%  2500  iterations) 
before  the  boundary  conditions  were  changed  to  initiate 
blowing.  The  same  anomalous  velocities  were  encountered 
near  the  slot. 

With  the  converged,  non-blowing  case  as  an  initial 
condition,  the  blowing  was  initiated  at  a  lower  momentum 
coefficient.  By  reducing  the  momentum  coefficient,  it  was 
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Figure  7.  Abnormal  velocities  at  slot  exit. 


Re=40 
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hoped  that  the  flow  out  the  slot  would  not  "overpower"  the 
freestream  flow,  Re=40,  M,=0.2.  Again,  the  velocity  vectors 
in  the  vicinity  of  the  slot  were  similar  to  those  presented 
in  Figure  7. 

In  the  hopes  of  solving  the  "overpowering"  problem  and 
to  investigate  more  realistic  conditions,  the  Reynolds 
number  was  increased.  Even  though  the  low  Reynolds  number 
regimes  are  not  practical,  they  were  to  be  compared  with 
potential  flow  theory  and  theoretical  procedures  developed 
for  jet  blowing  (8,27).  By  increasing  the  Reynolds  number, 
comparisons  with  experimental  data  can  be  made  (8,14,15)  as 
well  as  comparisons  with  potential  flow  theory  and 
theoretical  models.  However,  increasing  the  Reynolds  number 
introduces  many  complications.  In  the  subcritical  to  upper 
transition  regimes,  see  Appendix  D,  many  things  happen  to 
the  flow,  from  the  disappearance  of  a  regular  Karmen  vortex 
street  to  phenomenon  associated  with  boundary  layer 
transition.  A  Reynolds  number  of  3  x  10*  was  chosen  for  a 
variety  of  reasons:  to  avoid  the  boundary  layer 
instabilities  (i.e.  separation  bubbles)  associated  with  the 
subcritical,  critical,  and  supercritical  regimes;  to 
correspond  to  the  Reynolds  numbers  used  during  wind  tunnel 
testing  of  the  airfoil  (1);  and  to  correspond  to  the 
anticipated  wind  tunnel  testing  of  a  cylinder  at  AFIT. 

These  higher  Reynolds  numbers  required  the  inclusion  of  the 
Baldwin-Lomax  turbulence  model .  Based  on  experimental 
results  (2:631),  the  boundary  layer  transition  point  was 


36 


specified  as  sixty-five  degrees  from  the  front  stagnation 
point  for  both  upper  and  lower  surfaces.  The  blowing  was 
not  initiated  because  the  non-blowing  solution  never 
"converged"  to  a  steady-state  or  periodic/quasi-periodic 
state. 

To  stay  away  from  instabilities  caused  by  boundary 
layer  transition,  the  transition  point  was  specified  as  the 
front  stagnation  point.  Again  the  solution  never  converged. 
The  drag  coefficient  oscillated  between  0.7  and  0.9,  see 
Figure  8,  without  any  defined  period  (the  lift  and  moment 
coefficients  were  essentially  constant  at  a  value  of  zero). 
Figures  9  and  10  show  the  vortices  and  the  recirculation 
immediately  downstream  of  the  slot;  note,  the  slot  is 
located  between  0.49  and  0.50  on  the  vertical  axis.  The 
rear  separated  wake  region  is  no  longer  characterized  by 
symmetric  recirculation  vortices,  see  Figure  11. 

A  modification  to  the  Baldwin-Lomax  turbulence  model 
was  then  included  to  limit  the  distance  searched  for  maximum 
vorticity  in  the  rear  stagnation  region.  The  boundary  layer 
thickness  calculated  by  the  Baldwin-Lomax  model  on  the  front 
of  the  cylinder,  before  separation,  was  reviewed.  Based  on 
this  thickness,  a  limit  was  imposed  on  the  normal  distance 
the  Baldwin-Lomax  model  searched  for  the  maximum  vorticity 
(see  Appendix  C,  input  card  12).  This  modification  did  not 
create  any  major  changes  in  the  flow  field  results. 
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Figure  8.  Cylinder  Drag  Oscillations 


To  eliminate  any  additional  perturbations,  a  new  grid 
was  developed  for  a  cylinder  without  a  slot.  The  same  case 
was  run,  with  a  Reynolds  number  of  3  x  10s  and  the  boundary 
layer  transition  point  at  the  leading  edge.  The  same  basic 
results  were  obtained.  The  next  step  was  to  generate  a  new 
grid  with  more  points  in  the  boundary  layer  and  a  larger 
far-field  radius  (Awal 1= . 00001 ,  far-field  radius=20).  Again 
this  did  not  produce  a  "converged"  solution. 

A  variety  of  factors  may  contribute  to  the  lack  of  a 
steady  converged  solution  -  a)  actual  flow  phenomena,  b) 
turbulence  model,  and  c)  the  grid, 
a)  Flow  Phenomena 

The  Reynolds  number  examined  corresponds  to  the  upper 
transition  regime  (see  Appendix  D).  For  this  regime, 
experimental  tests  show  that  the  mean  lift  coefficient  is 
slightly  greater  than  zero  (0.02  to  0.03)  and  the  drag 
coefficient  is  increasing  from  approximately  0.2  to  0.5 
(22).  The  computed  lift  values  were  approximately  zero 
while  the  drag  coefficient  was  a  bit  high  at  0.8  to  0.9. 
However,  the  big  factor  is  the  lack  of  a  definitive  peak  in 
the  power  spectra  of  the  lift  fluctuations,  see  Figure  39. 
Schewe  (22)  plots  experimental  results  of  the  lift 
fluctuations  versus  the  Strouhal  number  (a  non-dimensional 
frequency)  for  a  variety  of  Reynolds  numbers.  For  the  upper 
transition  regime  of  Reynolds  numbers,  the  Strouhal  number 
is  characterized  by  broadband  low  frequency  values.  This 
broad  spectrum  implies  a  lack  of  periodicy  in  the  lift 
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fluctuations  and  the  vortex  shedding  frequency.  This  could 
be  responsible  for  the  computational  solutions'  erratic 
behavior.  Other  computational  investigations  have  cited 
this  as  a  reason  for  the  erratic  drag  oscillations  (10). 
b)  Turbulence  Model 

To  check  the  effect  of  the  turbulence  model,  iterations 
were  run  without  the  turbulence  model.  Starting  with  a 
uniform  flow  field,  the  code  was  run  for  100  iterations 
without  the  turbulence  model.  This  was  done  to  allow  the 
uniform  flow  field  to  adapt  to  the  cylinder  geometry.  At 
iteration  101,  the  Baldwin-Lomax  turbulence  model  was  turned 
on  (with  the  limit  on  the  search  for  maximum  vorticity  in 
the  normal  direction  and  without  the  Bradshaw  curvature 
modification).  After  approximately  4000  iterations,  the 
results  were  typical  of  previous  runs.  At  this  point,  the 
turbulence  model  was  tttrned  off.  After  approximately  3000 
more  iterations,  the  results  were  reexamined.  The  values 
for  the  drag  coefficient  were  not  as  erratic  or  as  large  in 
amplitude  as  the  values  obtained  with  the  turbulence  model. 
An  examination  of  the  parameters,  i.e.  FfI^e,  calculated  in 
the  turbulence  model  showed  a  rather  abrupt  change  upon 
reaching  the  wake  region. 

Obviously,  the  zero  equation  Baldwin-Lomax  turbulence 
model  is  not  the  best  choice  for  use  with  a  circular 
cylinder  with  large  separated  flow  regions.  Ideally,  the 
Baldwin-Lomax  model  is  most  easily  applied  to  a  C-grid 
configuration,  but  the  primary  problem  is  the  poor  handling 
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of  near  wake  regions  containing  separated  flow.  The 
circular  cylinder  precluded  the  use  of  a  C-grid  and 
amplified  turbulence  model  instabilities  due  to  the 
extensive  separated  wake  region.  Baldwin-Lomax  models  can 
not  accurately  model  regions  of  recirculating  flow.  For 
recirculating  flows.  Launder  and  Spalding  (13)  show  that  a 
two  equation,  eddy-viscosity  turbulence  model  (k-e  model)  is 
better  suited  for  these  regions  then  a  Baldwin-Lomax  model, 
c)  Grid 

For  high  Reynolds  numbers,  an  extensive  grid  is 
required  in  the  rear  quadrant  of  the  cylinder  to  properly 
model  the  flow.  Without  a  lot  of  grid  points  imbedded  in 
the  boundary  layer  and  near  wake  region,  all  of  the  flow 
idiosyncrasies  can  not  be  distinguished.  In  addition,  the 
far-field  location  is  critical  to  convergence.  Without  a 
significant  far-field  radius  the  wake  will  not  be  modeled 
correctly.  Another  major  consequence  of  an  inadequate 
far-field  location  is  the  reflection  of  pressure  waves. 

With  an  insufficient  far-field  location,  reflected  pressure 
waves  can  corrupt  the  solution. 

The  final  grid  used  for  this  case,  207  x  151,  is 
believed  to  have  been  fine  enough  to  resolve  the  above 
mentioned  problems.  The  normal  spacing  at  the  surface  was 
0.00001  and  the  far-field  boundary  was  circular  with  a 
diameter  of  30  chords.  The  same  results  obtained  with  the 
previous  grid  were  obtained  with  this  expanded  grid. 
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3.2  103RE  Airfoil 


As  stated  in  the  Introduction,  the  objective  was  to 
obtain  correlations  between  the  empirical  curvature 
correction  constant  and  the  flow  coefficients  (Cj,  C^,  and 
CB)  and  flow  field  around  a  circular  cylinder.  These 
correlations  were  then  to  be  compared  to  the  curvature 
correction’s  effect  on  the  flow  coefficients  and  flow  field 
around  the  airfoil.  However,  as  previously  covered,  a 
"converged"  solution  for  flow  around  the  cylinder  was  not 
obtained.  Instead,  effects  of  the  curvature  correction  on 
the  flow  around  the  streamlined  airfoil  were  investigated. 

In  previous  research,  a  range  of  values  for  the  Bradshaw 
curvature  correction  have  been  used.  Shrewsbury  obtained 
"results  more  nearly  in  agreement"  with  a  Bradshaw  value  of 
25  (24).  A  Bradshaw  value  of  3.2  "matched"  lift  coefficient 
values  but  not  drag  and  moment  coefficient  values  for 
Williams  (31).  After  varying  the  Bradshaw  constant  from  0  to 
10,  the  value  of  3.2  was  chosen  because  the  computed  lift 
coefficient  happened  to  coincide  with  the  experimental  lift 
coefficient  at  8=3.2,  for  a  single  momentum  coefficient. 

Since  the  drag  and  moment  coefficients  did  not  match  at 
8=3.2,  does  the  computed  flow  field  match  experimental 
observations?  Instead  of  just  "matching"  values,  can  a 
correlation  between  the  Bradshaw  curvature  correction  and  the 
computed  flow  field  be  made?  Using  the  same  data  points  used 
by  Williams,  see  Table  1,  the  Bradshaw  constant  was  varied 
from  0  to  25  to  monitor  the  effects  of  8  on  the  flow  field. 


45 


Table  1.  Computational  and  Experimental  Test  Conditions  (31) 


Pt 

M, 

Re  s 
(1x10® ) 

aeff 

P0/P« 

VT, 

CP^eiP 

35 

0.3 

3.11 

-0.92 

1.137 

0.956 

.0094 

36 

ED 

3.09 

-1.66 

1.284 

0.934 

.0179 

With  a  Bradshaw  value  of  zero.  Figures  12  and  13  show 
the  flow  field  characteristics  in  the  trailing  edge  region 
for  data  point  35.  The  Mach  contours.  Figure  12,  show  the 
freestream  flow  being  entrained  by  the  jet  and  separation 
occurring  just  prior  to  the  trailing  edge  point.  For  this 
low  of  a  blowing  coefficient,  Cp]c01p=O . 0094 ,  the  separation 
point  is  not  expected  to  reach  the  lower  surface  of  the 
airfoil.  Figure  13  shows  a  close-up  of  the  velocity  vectors 
in  the  separation  region.  The  separation  point  and  reverse 
flow  due  to  the  adverse  pressure  gradient  are  easily 
distinguishable.  Note  the  absence  of  any  vortex  in  the  flow 
field. 

As  the  Bradshaw  value  is  increased,  with  the  jet 
pressure  and  temperature  ratios  held  constant,  several 
important  flow  field  changes  occur.  The  separation  point 
tends  to  move  closer  to  the  jet  slot  and  two  counterrotating 
vortices  start  to  appear.  In  Figure  14,  with  a  Bradshaw 
value  of  4,  the  separation  point  has  moved  up  to  a  y/c 
location  of  approximately  0.010  versus  a  location  of  0.006 
in  Figure  13.  In  addition,  two  vortices  have  been  formed  in 
the  separation  region;  an  upper  vortex  with  a  clockwise 
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rotation  and  a  lower  vortex  with  a  counterclockwise  rotation 
against  the  airfoil  surface.  By  the  time  the  Bradshaw 
constant  is  10,  the  separation  point  is  above  0.012  and  two 
definite  vortices  exist  in  the  separation  region,  see 
Figure  15. 

With  8=25,  Figures  16  and  17  show  the  Mach  contours  and 
velocity  vectors,  respectively,  around  the  Coanda  surface. 

A  comparison  of  Figures  13  and  17,  show  the  substantial 
rearward  progression  of  the  separation  point  from 
approximately  y/c  =  0.006  to  y/c  >  0.012.  The  only  input 
parameter  changed  was  the  empirical  Bradshaw  curvature 
correction;  so  what  flow  pattern.  Figure  13  or  Figure  17,  is 
the  actual  model  of  the  physical  flow  field?  In  reference 
25,  Shrewsbury  diagrams  the  computational  and  experimental 
streamlines  for  a  higher  jet  momentum  value  and  for  a 
different  airfoil  geometry,  see  Figures  18  and  19.  However, 
these  figures  can  still  provide  insight  into  the  true  flow 
field  pattern.  Comparison  of  the  experimental  streamlines. 
Figure  19,  with  the  previously  discussed  computational 
streamlines.  Figures  12-17,  show  the  higher  Bradshaw  values 
provide  a  more  accurate  model  of  the  flow  field.  However, 
the  definitive  stagnation  streamline  separating  the  two 
vortices  is  not  evident  in  the  computational  solutions. 

A  better  model  of  the  flow  in  the  trailing  edge  region 
is  expected  if  a  more  sophisticated  turbulence  model  were 
employed.  Launder  and  Spalding  (13)  compare  a  two  equation 
eddy  viscosity  model  with  other  turbulence  models  for  a 
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Figure  15.  Velocity  vectors,  data  point  35,  Bradshaw=10 


Figure  12. 


Shrewsbury's  computational  streamlines 


Figure  13. 


Shrewsbury's  experimental  streamlines 


variety  of  flows.  Two  of  the  more  interesting  flow 
conditions  are  a  wall  jet  and  a  recirculation  zone.  Launder 
and  Spalding  show  the  k-e  model  provides  a  more  accurate 
prediction  of  the  velocity  decay  of  a  wall  jet  over  a  convex 
surface.  In  addition,  they  point  out  the  accurate  results 
obtained  with  the  two  equation  model  for  recirculation 
regions.  They  conclude  by  saying  the  two  equation  eddy 
viscosity  model  "is  the  simplest  kind  of  model  that  predicts 
both  near-wall  and  free -shear-flow  phenomena ... its  use  has 
led  to  accurate  predictions  of  flows  with  recirculation  as 
well  as  those  of  the  boundary  layer  kind."  For  these 
reasons,  the  adaptation  of  a  k-€  model  would  most  probably 
provide  better  predictions  of  the  flow  field. 

The  variation  of  the  Bradshaw  constant  has  an  effect  on 
more  than  just  the  flow  characteristics.  With  increasing 
Bradshaw  values,  the  various  flow  coefficients,  Cj,  ,  CB, 
and  Cp  change.  The  change  in  the  coefficients  with  Bradshaw 
constant  are  plotted  in  Figures  20-23.  The  experimental 
values  are  those  obtained  by  Abramson  for  the  103RE  airfoil 
(1).  All  of  the  coefficients  decrease  in  magnitude  as 
Bradshaw  is  increased.  But  the  more  important  trend  is  the 
asymptotic  behavior  of  the  curves.  Above  a  Bradshaw  value 
of  11,  all  of  the  parameters  tend  to  stabilize  at  some 
constant  value.  This  corresponds  to  the  drastic  flow  field 
change  between  8=0  and  10  and  the  relative  steadiness  of  the 
flow  field  between  8=10  and  25. 


55 


1.137 


56 


|  t  i  i  i  |>-i  i  rrpTT  i  |  i  i  i  i  ]  i  i  i  i  |  l  m  !  |  i  !  m  |  I  M  I  |  M  I  I  |  I  M  I  | 


cn 

00 

to 

m 

^4" 

m 

CM 

x— 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

6 

o 

d 

o 

o 

T3 

O 

o 

o 

o 

o 

o 

57 


Bradshaw  constant 

Figure  21.  Cd  versus  Bradshaw,  data  point  35 
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These  variations  in  the  flow  coefficients  bring  several 
questions  to  mind.  Could  more  accurate  flow  coefficients  be 
computed  if  the  momentum  coefficient  was  matched  with  the 
experimental  value  instead  of  matching  the  jet  pressure  and 
temperature  ratios?  What  value  of  jet  ratios  or  momentum 
coefficient  would  be  required  to  have  the  stable  region  of 
the  flow  coefficients  (0  >  11)  accurately  compare  with  the 
experimental  results?  Could  all  three  of  the  coefficients, 
lift,  drag,  and  moment,  aline  with  the  experimental  values 
for  the  same  jet  momentum  or  jet  pressure  ratios?  Figures 
24-26  answer  these  questions. 

On  Figures  24-26  are  four  lines  corresponding  to  the 
experimental  data  and  the  various  computational  variations. 
The  three  computational  lines  are:  1)  the  jet  pressure  and 
temperature  ratios  are  matched  with  experimental  data;  2) 
the  momentum  coefficient  is  matched  with  experimental  data; 
3)  a  constant  momentum  coefficient  that  best  predicts  the 
experimental  values  in  the  stable  region  (0  >  11).  As 
shown,  the  matching  of  jet  momentum  coefficients  provides  a 
better  prediction  of  the  actual  flow  values,  but  still  not 
perfect.  The  momentum  coefficient  had  to  be  increased  by  34 
percent,  to  C^O.0126,  before  computed  coefficients  came 
into  agreement  with  experimental  coefficients. 
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Bradshaw  constant 

Figure  26  Cm  versus  Bradshaw,  data  point  35 


The  Bradshaw  constant  did  not  have  a  noticeable  effect 


on  any  other  flow  characteristics.  The  boundary  layer 
transition  point  remained  constant  along  with  the  leading 
edge  stagnation  point.  Various  additional  plots  for  data 
point  35  are  contained  in  Appendix  E. 

The  next  step  is  to  see  if  these  same  trends  are 
evident  for  different  flow  conditions.  Going  to  the  higher 
jet  blowing  associated  with  data  point  36,  the  same 
variation  of  the  Bradshaw  constant  was  repeated. 

The  Mach  contours  and  velocity  vectors  for  0=0  are 
plotted  in  Figures  27  and  28.  As  is  expected  the  separation 
point  is  further  from  the  jet  slot  and  a  higher  degree  of 
entrainment  is  evident  when  compared  with  data  point  35. 

The  velocity  vector  diagram  shows  two  vortices  have  already 
been  formed.  Does  this  signify  a  breakdown  in  the  trends 
established  for  data  point  35?  Not  necessarily,  for  a 
Bradshaw  value  of  4,  see  Figures  29  and  30,  the  separation 
point  moves  toward  the  jet  slot  and  the  vortex  pattern 
changes  slightly.  The  y/c  location  of  the  separation  point 
for  0=4  is  located  at  approximately  0.005  compared  to  -0.005 
for  0=0.  As  seen  on  the  Mach  contour  plot,  the  strength  of 
the  upper  vortex  increases.  In  addition,  the  lower  vortex 
becomes  more  circular  in  nature  than  the  elongated  vortex 
for  a  Bradshaw  value  of  0.  These  same  general  trends 
continue  through  0=23,  see  Figures  31  and  32. 
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Figure  31.  Mach  contours,  data  point  36,  Bradshaw=23 
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Again,  the  question  arises  as  to  which  flow  pattern 
best  represents  the  real  flow  conditions?  Based  on  the 
previous  data  point  and  the  trends  with  this  data  point,  it 
is  safe  to  conclude  the  higher  Bradshaw  values  are  a  better 
indicator  of  the  true  flow  conditions. 

The  same  general  trends  for  the  flow  coefficients  are 
seen  in  Figures  33  -  36.  However,  the  Bradshaw  value  where 
the  coefficients  appear  to  reach  a  steady  value  is  increased 
to  15  or  more.  The  accuracy  of  the  lift  and  moment 
coefficients  is  comparable  to  the  lower  case,  but  the 
drag  coefficient  is  substantially  less  accurate.  Note,  the 
computed  Cj  coincides  with  experimental  results  at  0 
slightly  greater  than  two.  However,  the  flow  field  nor  the 
drag  and  moment  coefficients  agree  with  the  experimental 
results.  Appendix  F  contains  additional  plots  of  the  flow 
for  data  point  36. 
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versus  Bradshaw,  data  point  36 
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UD  Computational 
—  Experimental 


0.005 


IV.  COMCLUBIOHS 


Before  a  computational  algorithm  can  analyze  the 
benefits  of  a  concept,  it  must  first  be  able  to  accurately 
predict  the  performance.  This  report  did  not  try  to  analyze 
the  benefits  of  the  circulation  control  airfoil  concept. 
Instead,  the  purpose  of  this  study  was  to  investigated  the 
effects  of  the  empirical  Bradshaw  curvature  correction  on 
the  predictive  capabilities  of  a  Beam-Warming  approximate 
factorization  algorithim. 

4.1  Circular  Cylinder 

Using  the  circular  cylinder  to  get  a  handle  on  the 
Bradshaw  constant  is  a  worthwhile  investigation.  However, 
the  wake  dominated  characteristics  of  the  cylinder  flow 
field  introduce  many  complexities  not  characteristic  of  the 
airfoil  flow  field.  The  main  complication  is  the  unsteady 
flow  for  most  Reynolds  number  regimes.  Indeed,  this  study 
was  unable  to  achieve  a  steady  "converged”  solution  for  the 
flow  around  the  cylinder.  From  the  research  accomplished, 
the  lack  of  convergence  might  be  attributed  to  a  combination 
of  three  major  factors: 

1)  Actual  flow  phenomena 

2)  Inadequate  turbulence  model 

3)  Grid 
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4.2  103RE  Airfoil 


The  Bradshaw  constant  has  a  pronounced  effect  on  the 
flow  around  the  airfoil.  The  Bradshaw  value  influences  both 
the  flow  field  pattern  and  the  flow  coefficients,  Cj,  Cj, 
and  C,,.  The  development  of  the  flow  field  in  the  rear 
separation  region  approaches  agreement  with  experimental 
observations  as  8  is  increased.  However,  a  primary  and 
secondary  vortex  develop  instead  of  two  counter  rotating 
vortices  separated  by  a  stagnation  streamline.  The  flow 
coefficients  stabilize  as  the  Bradshaw  value  increases. 

This  stabilization  corresponds  to  a  stabilization  in  the 
computed  flow  pattern.  Beyond  this  stabilizing  point, 
further  increases  in  6  have  little  effect  on  the 
coefficients  or  the  flow  pattern. 

These  same  general  trends  apply  to  higher  blowing 
coefficients,  with  a  couple  of  exceptions.  The  variation  in 
the  flow  pattern  is  not  as  pronounced  and  the  flow 
coefficients  take  longer  to  level  off.  Instead  of  the 
development  of  two  vortices,  the  higher  blowing  case  saw  the 
refinement  of  the  two  vortices  with  increasing  6.  The  point 
of  stabilization  for  the  flow  coefficients  occurred  at  a 
higher  value  of  8.  Any  correlation  between  this  delay  and 
the  flow  refinement  is  not  totally  obvious  for  this  case. 

When  the  jet  ratios  are  matched  with  experimental 
values,  the  momentum  coefficients  computed  by  the  algorithm 
are,  in  general,  lower  than  the  experimental  momentum 
coefficients.  By  matching  the  momentum  coefficients  instead 
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of  the  jet  ratios,  the  computed  flow  coefficients  are 
slightly  better.  However,  the  values  are  still  not 
predictive  of  the  actual  experimental  results.  A 
substantial  increase,  more  than  30  percent,  in  the  momentum 
coefficient  was  required  before  computational  values  for 
lift  and  drag  were  in  close  agreement  with  the  experimental 
results.  The  moment  coefficient  tended  to  lose  accuracy  as 
the  momentum  coefficient  was  increased. 

Finally,  a  predictive  ability  was  never  achieved  with 
the  use  of  the  Bradshaw  modified  Baldwin-Lcmax  turbulence 
model.  A  higher  order  turbulence  model  must  be  incorporated 
before  any  attempt  can  be  made  at  predicting  the  flow  around 
a  circulation  control  airfoil. 
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V.  RECOMMEMDAT I QMS 


For  future  research,  the  following  recommendations  are 

made: 

1)  Apply  a  higher  order  turbulence  model.  A  two 
equation  eddy  viscosity  model  will  preclude  the  need  for  an 
empirical  curvature  constant  and  improve  the  accuracy  of  the 
computational  algorithim.  The  jet  velocity  decay  and  the 
recirculating  flow  can  be  modeled  more  accurately  with  a  k-e 
turbulence  model  (13). 

2)  Conduct  a  more  indepth  investigation  into  the  flow 
around  a  cylinder  at  high  Reynolds  numbers.  See  if  a 
"converged”  solution  can  be  obtained  for  a  Reynolds  number 
where  the  Strouhal  number  has  a  definitive  peak. 

3)  Conduct  an  experimental  investigation  concentrating 
on  the  flow  field  in  the  separation  region.  What  form  do 
the  streamlines  take  and  how  do  they  change  with  a  change  in 
the  flow  conditions?  An  analysis  of  the  turbulence  in  the 
separation  region  can  help  with  the  development  of  CFD 
turbulence  models  (17:3). 

4)  Conduct  a  theoretical  study  of  the  circulation 
control  flow  field.  Continue  Smith  and  Dunham's  (8,27)  work 
on  a  cylinder  with  a  jet  slot  and  try  to  expand  this 
mathematical  model  to  a  circulation  control  airfoil. 
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APPENDIX  A.  Non-dimensional  Variables 


In  general ,  the  governing  equations  are  non- 
dimensional  i  zed  for  convenience.  By  non-dimensional izing, 
characteristic  parameters  such  as  Mach  number,  Reynolds 
number,  and  Prandtl  number  appear  in  the  equations.  And  the 
flow  variables  are  normalized  as  a  result  of  non- 
dimensional  izing  .  To  non-dimensionalize,  a  reference 
length,  L,  must  be  established  along  with  reference  values 
for  velocity,  temperature,  density,  and  viscosity.  In  this 
case  the  chord  is  chosen  as  the  reference  length  and  the 
freestream  values  for  velocity,  temperature,  density,  and 
viscosity  are  used  (3:619).  The  non-dimensional  parameters 
are  denoted  as  ()*.  Note,  in  the  main  body  of  the  report 
the  superscripts  are  dropped. 
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Using  the  perfect  gas  law  and  substituting  in  for  p,  the 
non-dimensional  pressure  equation  becomes 


p.vm 


VJ 


p  *T* 
Y 


Thus,  the  non-dimensional  freestr'im  pressure  becomes 


The  dimensional  form  of  Southerland's  law  is  written  as 


ar2/2 
|1=  r+q 

When  non-dimensional i zed  Southerland's  law  becomes 


=  JL  = 

Pm 


qr,/a| 

’  r„  +  c2 

T+C2  J 

C^8 

—  ^  y*)  3/2 


r.  +  c. 


or 


where 


(T*) 


1  +  C2* 

r*  +  Cl 
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The  parameter  k/cf  shows  up  often  in  the  governing 
equations.  The  non-dimensional  version  is  obtained  from  the 
definition  of  the  Prandtl  number 


Pr=£ &  = 

k  cv  k/cv 

Solving  for  k/cf 

k  _ 

cv~  Pr  Pr 

It  the  chord  length  (reference  length)  is  one  then 


thus 


*  -  vp* 

cv  Re  Pr 
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APPENDIX  B .  Jacobian  Matrices 
The  Jacobian  matrices  specified  in  Eqs  64  and  65  are 
expanded  below.  Before  calculating  the  derivatives,  the 
vectors  are  rewritten  in  terms  of  the  conserved  variables 
(6:98) 


A  = 


60 
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$y*-VUc 
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uc- (y-2)ixu 
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0 
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(Y-i)ny 
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The  R  and  S  matrices  are  defined  as 


where 


dV, 

dtfl 


f  0  0  0  0  } 

^  -(bxu+b2v)  jbx  b2  0 

~p  -  (b2u  +  b2v)  b2  i?3  0 

k  2?41  i?42  ^43  ^44, 


jb4i«-  [jbf  +  2b2uv  +  b2v2  +i>4y  (y-l)  (u2*^2-^)  ] 
^4a-i>i“^4Y<Y-l)  u  +  b2v 
Jb43*Jb2u+jb3-jb4Y(y-l) 

^44*jb4Y<Y“1) 
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where 


dW% 


0  0  0 

.  -(cfiU  +  dgV)  d1  <4 
S-  — 

p  -( d2u  +  d3v )  <4  <4 

d41  d42  3*3 


d41*  -  [di +2d2uv’  +  d3v2  +  d4Y  (Y”D  ( u2+v2-Et )  ] 
d^di-drfiy-Du+dzV 

d43  =  d2u+d3~d4Y(Y-D 
d44*d4Y(Y-D 
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The  Beam-Warming  code  developed  by  Dr  Miguel  Visbal 
allows  a  wide  variety  of  user  inputs.  The  user  can  specify 


the  type  of  scheme,  type  of  grid,  turbulence  model 
parameters,  time  stepping  'Employed,  jet  slot  conditions, 
among  others.  The  following  list  describes  all  of  the  user 
defined  inputs  (3i). 

CARD  1: 

H,  PHI  -  Time  differencing  scheme  to  be  employed 

with  the  Beam-Warming  algorithm.  H  and 
PHI  correspond  to  8j  and  02,  respectively, 
as  seen  in  most  documentation. 


Scheme 

H 

PHI 

Order 

Euler  explicit 

0 

0 

first 

Euler  implicit 

1 

0 

first 

Trapezoidal 

0.5 

0 

second 

3-point  backward 

1 

0.5 

second 

Leap  frog 

0 

0.5 

first 

CARD  2: 

NSTEPI  -  Iteration  number  to  start  with. 

TAU  -  Physical  time  at  NSTEPI. 


CARD  3: 

IL  -  Number  of  grid  lines  in  azimuthal  (or  wrap 

around)  direction.  If  any  grid  lines  are 
overlapped  make  sure  they  are  included  in 
the  value  for  IL.  If  using  the  INGFDLNS 
program  then  IL  is  four  more  than  the 
number  of  grid  lines  originally  built  due 
to  overlapping. 

JL  -  Number  of  grid  points  in  normal  direction. 

ILE  -  Grid  index  for  airfoil  leading  edge.  If 

INGPDLNS  has  been  used  than  ILE  will  be  3. 
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I  TEL 


Grid  index  for  the  airfoil  trailing  edge. 


ITEU 

CARD  4: 

XM1 

RE 

TW 

SI 

ALFA 

CARD  5: 

IADW 

CARDS  6-7: 

OMEGAO , 
TOMGO , 

AFLO , ALF1 . 
RFREQ 

CARD  8: 

CHOM 

XC,  YC 


This  parameter  is  for  an  upper  trailing 
edge  location.  However,  the  program  does 
not  use  this  parameter,  so  ITEU  =  0. 


Freestream  Mach  number. 

Reynolds  number  based  on  airfoil  chord. 

-  Hall  temperature  normalized  by  freestream 
temperature  (see  IADH  parameter). 

Parameter  for  Southerland's  law  (SlsCj*). 

-  Angle  of  attack  (degrees). 


Flag  for  adiabatic  wall  conditions.  If 
IADH  is  0  then  TW  must  be  specified;  if 
IADH  is  1  then  TH  must  be  999. 


These  parameters  are  for  non-stationary 
airfoils.  For  analysis  with  stationary 
airfoils  (as  done  in  this  study)  these 
parameters  are  0. 


-  Reference  point  in  the  x-direction, 
normalized  by  airfoil  chord,  for 
computation  of  the  pitching  moment. 

-  These  parameters  are  for  non-stationary 
airfoils.  For  analysis  with  stationary 
airfoils  (as  done  in  this  study)  these 
parameters  are  0. 
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CARD  9: 


I0FLW1 


IOFLW2 

CARD  10: 

IMOVE, 

IMPB 

CARD  11: 
NTURB 


I  FREQ 


CARD  12: 
JLAST 


Inflow/outflow  boundary  transition  point 
for  the  upper  surface.  From  this 
azimuthal  grid  index  to  IOFLW2,  the 
outflow  boundary  conditions  are  imposed  at 
the  far-field  boundary.  A  typical  value 
is  a  grid  index  3-5  points  past  midchord. 

Inflow/outflow  boundary  transition  point 
for  the  lower  surface. 


These  parameters  are  for  non-stationary 
airfoils.  For  analysis  with  stationary 
airfoils  (as  done  in  this  study)  these 
parameters  are  0. 


Iteration  number  for  inclusion  of  the 
turbulence  model.  For  high  Reynolds 
numbers,  the  turbulence  model  should  be 
included  when  DTVIS  is  greater  than  zero. 
Typically,  The  turbulence  model  is 
included  in  the  algorithim  at  iteration 
101. 

Once  the  turbulence  model  is  included  in 
the  analysis,  I FREQ  specifies  now  often 
the  turbulence  is  recalculated. 

Nominally,  the  turbulence  model  should  be 
recalculated  every  iteration,  so  IFREQ  = 
1. 


This  is  a  modification  to  the  Baldwin- 
Lomax  turbulence  model  to  handle  wake 
regions.  It  is  a  limit  on  how  far  out 
(normal  direction)  the  turbulence  model 
will  search  for  the  end  of  the  boundary 
layer  region.  JLAST  should  be  set  to  a 
value  corresponding  to  the  JMAX  values 
calculated  in  the  flow  field  before 
separation. 
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JFMAX  -  This  is  associated  with  the  previous 

parameter.  This  parameter  limits  the 
number  of  indices  in  the  normal  direction 
to  find  FMAX  and  JMAX.  Nominally,  JFMAX  = 
JLAST  -  1.  In  addition,  the  output  flag 
I FLAG  must  be  watched.  If  IFLAG=0 ,  the 
turbulence  model  is  modeling  a  normal 
boundary  layer,  otherwise  the  turbulence 
model  is  modeling  a  separated  wake  region. 

IJ1,  IJ2  -  Azimuthal  grid  indices  for  a  region  (on 

the  surface)  where  the  turbulence  model  is 
not  used.  This  option  is  intended  for 
use  around  sharp  edges,  i.e.  jet  slot. 

I JBL  -  This  defines  the  type  of  boundary  layer 

for  the  region  between  IJ1  and  IJ2.  If 
I JBL  is  set  to  0,  then  the  region  is 
laminar.  If  IJBL=1,  then  the  region  is 
turbulent  and  the  eddy  viscosity  is 
linearly  interpolated. 


CARD  13: 

ITRAN  -  Iteration  number  when  the  code  starts  to 

use  the  pressure  gradient  transition 
criteria.  With  this  criterion  the 
boundary  layer  transition  point  will  be  at 
the  point  of  minimum  pressure,  if  the 
point  of  minimum  pressure  is  located 
between  the  stagnation  point  and  DSTR1/2. 
This  criteria  will  give  a  conservatively 
high  estimate  of  the  drag.  This  criteria 
should  not  be  used  until  the  turbulent 
boundary  layer  has  converged  sufficiently. 

DSTR1  -  Upper  surface  boundary  layer  transition 

location.  Measured  in  arc  lengths, 
normalized  by  the  chord,  form  the 
stagnation  point.  See  reference  23  Fig 
17.9  for  typical  values. 

DSTR2  -  Lower  surface  boundary  layer  transition 

location. 
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CARD  14: 


NDTAU  -  This  number  specifies  the  iteration  when  a 

new  time  step  will  be  computed.  Typically 
this  option  is  not  exercised. 

CFL  -  Courant-Fredrichs-Lewy  number.  To  avoid 

instabilities,  CFL  should  be  1  for  the 
first  several  iterations  (approximately 
100)  while  the  initial  condition  is 
adapting  to  the  current  geometry.  Later 
CFL  can  be  increased,  typically  to  10. 

BETA  -  Specifies  the  type  of  time  stepping  used 

by  the  algorithim. 

BETA=0  --  time  accurate  solution  where 
the  time  step  is  constant 
and  equal  to  DTVIS. 

BETA=1  —  local  time  stepping  for 
steady  state  solution. 

The  time  step  is  the  maximum 
of  DTVIS  or  CFL  times  DTCFL . 

DTVIS  -  Time  interval.  For  local  time  stepping 

DTVIS  should  be  0  for  initial  iterations. 
For  later  iterations  set  DTVIS  to  .002  if 
possible. 


CARD  15: 

ISPECT  -  If  set  to  one  a  scaling  factor  will  be 

calculated  for  the  fourth  order  damping 
coefficients.  If  ISPECT=0  then  no  fourth 
order  damping  will  be  added.  Typically 
ISPECT=1. 

WE,  WI  -  Fourth  order  damping  factors.  Nominally 

WE  should  be  between  .01  and  .02  with  WI 
between  .025  and  .05.  Typically  WE=.01 
and  WI=.025.  If  these  parameters  have  a 
significant  effect  on  the  solution  then 
the  grid  resolution  may  not  be  sufficient 
for  the  case  being  examined. 

WPE,  WPI  -  Additional  damping  factors  to  be  used  in 

regions  where  supersonic  flow  exists. 
Should  be  between  .1  and  .5.  For  flows 
without  any  supersonic  regions  these 
values  should  be  0. 


90 


CARD  16: 


INMAX 

NCONV 

CARD  17: 

ISLOT1 , 
ISLOT2 

TPR 


TTR 


CARD  18: 
ICURV 

CALP 

FCMIN , 
FCMAX 


Maximum  number  of  iterations  to  perform. 

Convergence  and  force  data  is  written  to 
appropriate  files  every  NCONV  iterations. 


Azimuthal  indices  defining  the  location 
of  the  jet  slot. 

The  non-dimensional  ratio  of  jet  total 
pressure  to  freestream  static  pressure. 


The  ratio  of  jet  total  temperature  to 
freestream  static  temperature. 


If  this  flag  is  0,  then  the  Bradshaw 
curvature  correction  is  not  used.  If 
ICURV  is  1,  then  the  Bradshaw  curvature 
correction  is  used. 

Is  the  empirical  constant,  0,  in  the 
Bradshaw  curvature  correction  equation. 

-  The  minimum  and  maximum  values  of  the 

Bradshaw  curvature  correction.  Nominally, 
values  of  0.5  and  1.5,  respectively,  are 
used. 
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APPEHDIX  D.  Plow  Regimes  Around  a  Cylinder 


Even  though  the  geometry  may  be  "simple",  the  flow 
around  a  two-dimensional  circular  cylinder  is  very 
complicated.  The  following  list  is  a  breakdown  of  the 
various  flow  regimes  around  a  cylinder.  The  information  was 
obtained  form  a  variety  01.  sources  and  the  various  regime 
titles  are  those  most  commonly  used  throughout  the 
literature  (2,21,22).  The  Reynolds  number  breakpoints  are 
approximate;  exact  values  are  a  function  of  freestream 
turbulence,  surface  roughness,  etc. 


Reynolds  number 
0-5 


5-40 


40  -  80 


80  -  5000 


Flow  Description 

"Stokes  range".  This  regime  is  typified 
by  no  separation  and  a  seemingly 
inviscid  fluid.  The  drag  coefficient  is 
very  high  but  drops  off  quickly  as  the 
Reynolds  number  is  increased,  see 
Figure  37. 

"Symmetrical  wake  region".  The  flow 
separates  from  the  cylinder  with  the 
appearance  of  two  counterrotating 
vortices.  The  flow  is  laminar  up  to  the 
separation  point.  The  pair  of  captive 
vortices  are  symmetric  about  the 
cylinder  midplane.  As  the  Reynolds 
number  is  increased  the  vortices 
elongate.  The  drag  coefficient  is 
steadily  decreasing. 

"Incipient  Karmen  vortex  range". 
Instability  sets  in  causing  the  vortices 
to  become  asymmetric  and  ultimately 
shedding  occurs. 

"Pure  Karmen  range".  A  regular  Karmen 
vortex  street  is  established.  The 
shedding  frequency  becomes  periodic  and 
is  characterized  by  the  Strouhal  number. 
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5000  -  2  x  10*  "Subcritical  range".  The  disappearance 

of  a  regular  vortex  street  and  more  of  a 
quasi-periodic  shedding  typify  this 
regime;  the  Strouhal  number  is  constant 
at  approxim  rely  0.2.  Laminar 
separation  occurs  before  the  flow 
reaches  the  vertical  mid-plane. 

2  -  4  x  10*  "Critical  range".  In  this  regime,  the 

drag  is  sharply  reduced  due  to  boundary 
layer  transition.  The  laminar  boundary 
layer  separates  from  the  cylinder. 
Immediately  after  separation,  transition 
occurs.  The  turbulent  flow  reattaches 
to  the  body  leaving  a  separation  bubble 
on  the  surface  of  the  cylinder.  The 
higher  energy  of  the  turbulent  boundary 
layer  continues  to  flow  along  the  body 
before  finally  separating  near  the  aft 
end.  The  drag  is  reduced  due  to  the 
decrease  in  the  wake  region.  At  first 
the  flow  is  asymmetric  due  to  the 
existence  of  only  one  separation  bubble. 
As  the  Reynolds  number  is  increased,  a 
separation  bubble  forms  on  the  opposite 
side  and  flow  symmetry  returns.  The 
vortex  shedding  frequency  is  still 
quasi-periodic.  The  Strouhal  number  is 
no  longer  constant  but  has  sharp 
discontinuities  brought  on  by  the  flow 
phenomenon,  see  Figure  38. 

4  -  10  x  105  "Supercritical  range”.  The  separation 

bubbles  stabilize  along  with  the  drag 
coefficient.  The  quasi-periodic 
shedding  continues. 

1  -  5  x  10*  "Upper  transition".  The  quasi-periodic 

sheddir g  is  replaced  by  unperiodic 
broadband  low  frequency  shedding,  see 
Figure  39.  The  boundary  layer 
transition  occurs  without  the  presence 
of  the  separation  bubbles  and  the  drag 
coefficient  rises  from  the  supercritical 
plateau  (0.2)  to  the  transcritical 
plateau  (0.5). 

5  x  10*  and  up  "Transcritical  range".  The  periodic 

vortex  shedding  reappears  in  the 
transcritical  regime  and  the  drag 
coefficient  levels  off  at  0.5. 
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Figure  37  represent  effect  of  Reynolds  number  on  the 
drag  coefficient.  Figure  38  is  a  closer  look  at  the 
supercritical  to  transcritical  regimes.  The  power  spectra 
of  the  lift  fluctuations  is  outlined  in  Figure  39.  Note  the 
lack  a  defined  peak  in  pictures  c  and  d. 
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Figure  37.  versus  Reynolds  number  for  a  cylinder  (23) 


Figure  39.  Power  spectra  of  Lift  Fluctuations, 
(b)  Re=7.2xl0'  (c)  Re=l. 9x10s,  (d)  Re=3.7xl0\  ( 
(f)  Re=7 . 1x10°  (ref  22). 


Appendix  E.  Data  Point  35  Plots 


98 


100 


Figure  35.  Velocity  vectors,  data  point  35,  Bradshaw=8 


Figure  36.  Velocity  vectors,  data  point  35,  Bradshaw=14 
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Figure  37.  Velocity  vectors,  data  point  35,  Bradshaw=20 


Bradshaw=20,  Cl  =0.4369,  Cd=0. 007185,  Cmu=0. 00920 


Figure  38.  Velocity  vectors,  data  point  35,  Bradshaw=20 


Appendix  F. 
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Data  Point  36  Plots 
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mach  plt  •*  Bradshaw=2,  _CI  =  1 .126.  Cd=0  01 150,  Cmu=0  01834 


105 


Figure  39.  Mach  contours,  data  point  36,  Bradshaw=2 


Bradshaw=10,  CU0.8992,  Cd=0. 007701,  Cmu=0. 01830 


Figure  42.  Velocity  vectors,  data  point  36,  Bradshaw=10 


Bradshaw=14,  CU0.8740,  Cd=0. 007290,  Cmu 
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Figure  43.  Velocity  vectors,  data  point  36,  Bradshaw=14 


vel.plt  **  Bradshaw=20,  Cl=0. 8576,  Cd=0. 007067,  Cmu=0. 01844 


Figure  45.  Velocity  vectors,  data  point  36,  Bradshaw=20 
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